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■ Abstract 

o : 

I The improvement of simulations of QCD with dynamical Wilson 

' fermions by combining the Hybrid Monte Carlo algorithm with parallel 

. tempering is studied on lO'^ and 12^ lattices. As an indicator for 

I decorrelation the topological charge is used. 

CD . 

1 Introduction 

^> : 

^ i Decorrelation of the topological charge in Hybrid Monte Carlo (HMC) sim- 

^ [ ulations of QCD with dynamical fermions is a long standing problem. For 

■ - - ' staggered fermions an insufficient tunneling rate of the topological charge Q 

has been observed [jl|, ^. For Wilson fermions the tunneling rate has been 
reported to be adequate in many cases ||, |[. However, since the compari- 
son is somewhat subtle, the reason for this could also be that one is not as 
far in the critical region as with staggered fermions. In any case for Wilson 
fermions on large lattices and for large values of k near the chiral limit the 
distribution of Q is not symmetric even after more than 3000 trajectories 
(see e.g. Figure 1 of 0]). 

While many observables appear not sensitive to topology there are clearly 
ones for which proper account of topological sectors is important. In fact, 
it has been observed that the r]' correlator is definitely Q dependent 0. 
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Thus it is quite important to look for simulation methods that produce good 
distributions of Q. 

Generally one wants to penetrate as deeply as possible into the critical 
region. This, however, is limited by the increasing autocorrelation times. 
Thus better decorrelation is also desirable with respect to this goal. The 
topological charge appears to be a good touchstone of the improvements 
achieved by new methods. 

The method of simulated tempering has been proposed in where the 
inverse temperature has been made a dynamical variable in the simulations. 
The principle is, however, much more general in the sense that any parameter 
in the action can be made dynamical. The mechanism behind the better 
decorrelation of the simulation process is that in place of suppressed tunneling 
an easy detour through parameter space with little suppression is opened. 

In fact, considerable improvements have been obtained by simulated tem- 
pering with a dynamical number of the degrees of freedom in the Potts-Model 
0, with a dynamical inverse temperature for spin glass |§] and with a dy- 
namical monopole coupling in U(l) lattice theory |p. An investigation with 
a dynamical mass of staggered fermions in full QCD UTO) has indicated a 
potential gain at smaller masses. Simulated tempering requires the deter- 
mination of a weight function in the generalized action. Thus developing 
efficient methods |]^, for obtaining the weight function has been a crucial 
issue. 

A major progess was the proposal of the parallel tempering method |TT 



in which no weight function needs to be determined. This has allowed large 
improvements in the case of spin glass fill. In QCD also improvements 



have been obtained with staggered fermions |T2[. On the other hand, in 
simulations of QCD with 0(a)-improved Wilson fermions ||l3l no gain has 
been found. However, the use of only two ensembles in this study did not 
take advantage of the main idea of the method. 

In the present paper parallel tempering is used in conjunction with HMC 
to simulate QCD with (standard) Wilson fermions. To study the performance 
of the tempering method we have recorded the time series of the topologi- 
cal charge. In previous work we have done this on an 8^ lattice where 



already a considerable increase of the transitions of the topological charge 
was observed. However, the effect of enlarging the lattice size remained a 
major question. Furthermore, the rather narrow distribution of the topo- 
logical charge for this lattice size prevented us from resolving more details. 
Therefore, investigations on 10^ and 12^ lattices are done here. 
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2 Parallel tempering 



In standard Monte Carlo simulations one deals with one parameter set A and 
generates a sequence of configurations C . The set A includes /5, k and as 
technical parameters the leapfrog time step and the number of time steps 
per trajectory. C comprises the gauge field and the pseudo fermion field. 



In the parallel tempering approach [|TT| , |T5[ one simulates N ensembles 
(Aj; Cj), i = 1, . . . , in a single combined run. The whole run represents 
one bigger ensemble in an enlarged configurations space stratified (with re- 
spect to A) into the ensembles above. Two steps alternate: (a) update of 
configurations in the standard way, (b) exchange of configurations by swap- 
ping pairs. Swapping of a pair of configurations means 

((\ n\ (\ n w J ^i)' (-^ji ^«))' if accepted , ^ 

((A., C), (A„ C,)) j ^^^^ ^^^^^ ^^^^ (Ij 

with the Metropolis acceptance condition 

Pswap(^,j)=min(l,e-^^) , (2) 

AH = H,^ (Q) + H,^ (C,) - H,^ (C,) - H,^ (Q) (3) 
where H is the Hamiltonian of the HMC dynamics. Since after swapping due 



to detailed balance both ensembles remain in any case in equilibrium |jTT], |T5 
the swapping sequence can be freely chosen. In order to achieve a high swap 
acceptance rate one will only try to swap (/?, K)-pairs that are close together. 
If the chosen {/3, K)-values lie on a curve in the (/?, K)-plane there are three 
obvious choices for the swapping sequence of neighboring [jS, K)-pairs. One 
can step through the curve in either direction or swap randomly. We have 
observed that it is advantageous to step along such a curve in the direction 
from high to low tunneling rates of Q. 



3 Simulation details 

Wilson fermions and the standard one-plaquette action for the gauge fields 
have been used. The HMC program applied the standard conjugate gradient 
inverter with even/odd preconditioning. The trajectory length was always 1. 
The time steps were adjusted to get acceptance rates of about 70%. In all 
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cases 1000 trajectories were generated (plus 50-100 trajectories for thermal- 
ization) . 

Q was measured using its naively discretized plaquette form after doing 
50 cooling steps of Cabibbo-Marinari type. This method gives close to integer 
values which were rounded to the nearest integers. 

Statistical errors were obtained by binning, i.e., the values given are the 
maximal errors calculated after blocking the data into bins of sizes 10, 20, 50 
and 100. 

4 Results 

We have performed tempered runs using 6 and 7 ensembles, all at /5 = 5.6 
on 10^ and 12^ lattices, as well as standard HMC runs at fixed n values 
for comparison. Our ensembles cover the K-range investigated by SESAM 
(«; = 0.156, 01565, 0.157, 0.1575) In the run using 6 ensembles we studied 
the effect of extending the K-range by adding lower values of k, while in the 
run using 7 ensembles we have tested the efficiency of a denser spacing of the 
K- values. Our K-values are listed in Table |1|. 

Figures 1 and 2 show time series of Q obtained on 10^ and 12^ lattices, 
respectively, for standard HMC and for tempered HMC with 6 and with 7 
ensembles. One sees that tempering makes Q fluctuating much stronger. 
Thus correlations between subsequent trajectories indeed decrease. Com- 
paring the time series of Q on lO'* and 12^ lattices, the increasing width 
of the topological-charge distribution is seen to lead to a richer pattern of 
transitions. 

On the 10^ lattice, in addition to the tempered run with 7 ensembles 
with finer steps in k and the one with 6 ensembles which penetrates deeper 
into the uncritical region we also performed a tempered simulation which 
used only the four common K-points of the latter ones. The results of this 
indicate that both of these modifications cause only moderate improvements, 
which suggests that optimizing range and distances one might do with fewer 
ensembles. 

In principle an autocorrelation analysis would allow to quantitatively as- 
sess the improvement of decorrelating. However, such an analysis cannot 
be carried out with the given size of samples. Some quantitative account is 
never less possible considering the mean absolute change of Q, called mobility 
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Ill 



^Q = ^E (4) 
and the HMC time between topological events Tq given by 

1/Tg = -— 5:min(l,|Q(^)-Q(^-l)|) . (5) 

While Mq measures the change in topological charge, Tq tells how often 
changes (without regard to their magnitude) occur. The difference between 
Mq and 1/Tq increases with the width of the distributions of the topological 
charge. 

The results obtained for Mq and 1/Tq are given in Table |T]. Comparing 
Mq for standard and tempered runs, gains by factors 2 to 5 are obvious. At 
larger k, values the improvement in the tempered run with 7 ensembles is 
larger than in that with 6 ensembles. Thus using finer steps in k pays off 
more than including additional K-values in the uncritical region. Comparing 
Mq and 1/Tq a significant effect of the broader Q-distribution can be seen. 
As expected, this effect becomes stronger on the larger lattice. Generally one 
observes that both Mq and 1/Tq increase with the lattice size. 

The results obtained for (Q) and for (Q^) are given in Table |^. The 
values for (Q) are in general better localized at zero in case of the tempered 
runs. The results for (Q^) increase with lattice size, as is to be expected, and 
also for smaller k. The errors of {Q"^) are rather large for the standard run 
and for the tempered run with 6 ensembles while for the tempered one with 
7 ensembles they are moderate. This again indicates that under the given 
conditions smaller spacing of k is more advantageous. 

With respect to these errors one has to realize that within a tempering run 
the values found for (Q^) at different k, are correlated. This phenomenon is 
well known from autocorrelation curves and mass determinations. To account 
for it in the fit to a corresponding function the full weight matrix is to be 
used. To obtain this matrix fortunately one can rely on the fact that higher 
than two-point cumulants vanish and thus express four-point correlations by 



their disconnected parts |16|. 



Since we do not intend to make a fit here it suffices to keep in mind that, 
while the errors account properly for the individual values, they do not tell 
about relative errors of neighboring values. In any case, better data from 
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tempering, though correlated in k, are to be preferred to poor ones from 
standard runs. 

Concerning the magnitudes of (Q^) in Table though errors are large, 
for the 12^ lattice some possible tendencies may be discussed. For large k 
the results of the tempered runs with 6 and with 7 ensembles agree and 
have larger values than the standard run. Since the latter, as is obvious 
from Figure 2, only occasionally escapes from the value Q = this can be 
understood. On the other hand, for smaller k the values obtained in the run 
with 7 ensembles are below those of the other runs. In that case in Figure 
2 the curve of standard HMC is seen to travel also to larger \Q\, however, 
also to take quite some time to get back. A broadening of the distribution 
caused by this behavior explains the deviation. Because of less accuracy and 
decorrelation, the run with 6 ensembles in its behavior appears to be closer 
to the standard HMC case. 



5 Swap acceptance rates 

With regard to large scale simulations of QCD performance predictions are 
needed. One potential problem of the tempering method has been addressed 
T3[] , namely the decrease of the swap acceptance rate (A) with the lattice 
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volume. There the relation 



{A)=eTk{-^{AH)) (6) 



has been used which has been shown in to hold for Metropolis-like al- 



gorithms. In |13[, in simulations with 2 ensembles, the relation (y) has been 
found to hold for a large range of (AH). 

Our swap acceptance rates (A) for the 10^ and 12^ lattices here and for 
the 8^ lattice in ill are given in Table 0. To check the relation (^) we have 
determined y in 

{A) = erfc(y) (7) 

and have also listed the resulting quantities lOOy/L^ and lOy/L for our 
lattices in Table ^ (the factors of 100 and 10 where introduced to obtain 
numbers of 0(1)). For linear scaling of (AH) with lattice volume by (^ one 
expects the quantity lOOy/L^ to be constant. From Table ^ it is seen that this 
only holds for the case of 6 ensembles; for the case with 7 ensembles, instead, 
the quantity lOy/L gets constant. In this context it should be remembered 
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that the data from 7 ensembles, i.e. from smaller Ak, are our more accurate 
ones. 

With respect to the data on the 8^ some caution is necessary because in 
that case, a.t (3 = 5.6 and 0.15 < k < 0.16, the finite temperature phase 
transition |jl8| is crossed. This could influence the scaling behavior of (AH). 
However, in Table |^ the 8^ data (within errors) fit into the trend of the 
other data. Thus, the small-volume phase transition does not affect our 
observation. 

Clearly more ensembles will be needed on larger lattices if one wants to 
keep (A) and the parameter range constant. However, firstly, as our example 
with 7 ensembles shows, for smaller Ak, one does better than is expected from 
naive scaling arguments. Secondly, an open question is how the decrease of 
(A) and the slowing down of tunneling between topological sectors compete. 
Within this respect one can hope that the speed-up more than compensates 
the additional effort by taking more ensembles. Thirdly, since results at 
several parameter values are needed anyway, having more points with less 
statistics is essentially equivalent to fewer points with more statistics. 



6 Discussion 

We have observed that parallel tempering considerably enhances the tunnel- 
ing between different sectors of the topological charge. This enhancement 
also indicates an improvement of decorrelation for other observables. The 
method is particularly economical when several parameter values have to be 
studied in order to analyse the physical dependence. These features make 
parallel tempering an attractive method for QCD simulations. 

So far we have used equidistant parameter values. As known from simu- 
lated tempering ^] and from parallel tempering for spin glass |Tl| refine- 



ments using optimized distances are possible. In the present case such opti- 
mizations have been not practicable because of the amount of computer time 
which this would have needed. Similarly optimizing the parameter range, 
i.e. the penetration into the range with easy tunneling, has not been feasible. 
The economics of computer time enforces to defer such optimizations to the 
applications, i.e., "to learn on the job". 

In large-scale QCD simulations the performance on larger lattices is im- 
portant. It has been seen here that definite predictions for this from smaller 
lattices are difficult. Nevertheless, for example, to make progress with the 
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vl problem, one should just try parallel tempering. Where data at several 
parameter values are needed, nothing can be lost doing the respective simu- 
lations simultaneously and making the parameters dynamical. 

On the other hand, there are, of course, other problems in QCD for 
which the method is attractive. Thus one might think of finite temperature 
studies, of phase structures with modified actions, or generally of questions 
where slowing down so far prevents from satisfactory solutions. 
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Table 1: Values of Mq and 1/Tq at /3 = 5.6. 



standard HMC tempered HMC 

6 ensembles 7 ensembles 









Ak 


— 0005 


Ak — 


00025 




10^ 




10^ 


12^ 


10^ 


12^ 


K 








71 /T 






U.iOOUU 






U.bUl^D ) 


0.1 (\p) 






0.15550 






0.64(5) 


0.90(8) 






0.15600 


0.09(2) 


0.30(4) 


0.48(5) 


0.80(8) 


0.51(4) 


0.91(5) 


0.15625 










0.57(4) 


1.20(7) 


0.15650 


0.07(2) 


0.33(3) 


0.36(5) 


0.75(7) 


0.48(5) 


1.17(7) 


0.15675 










0.38(4) 


1.12(10) 


0.15700 


0.08(2) 


0.20(5) 


0.31(5) 


0.53(7) 


0.33(5) 


0.96(9) 


0.15725 










0.23(5) 


0.82(7) 


0.15750 


0.05(1) 


0.09(2) 


0.15(3) 


0.28(4) 


0.14(4) 


0.51(7) 


n 














0.15500 


0.24(3) 


0.38(2) 


0.35(2) 


0.48(2) 






0.15550 






0.46(3) 


0.53(2) 






0.15600 


0.09(2) 


0.27(3) 


0.40(4) 


0.45(3) 


0.42(3) 


0.55(2) 


0.15625 










0.48(3) 


0.69(3) 


0.15650 


0.07(2) 


0.28(2) 


0.32(4) 


0.47(2) 


0.41(3) 


0.69(2) 


0.15675 










0.32(3) 


0.64(3) 


0.15700 


0.08(2) 


0.18(4) 


0.28(4) 


0.36(3) 


0.27(4) 


0.61(3) 


0.15725 










0.21(4) 


0.58(3) 


0.15750 


0.05(1) 


0.09(2) 


0.14(3) 


0.22(3) 


0.13(4) 


0.36(4) 
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Table 2: Values of {Q) and (g^) at /3 = 5.6. 



standard HMC tempered HMC 

6 ensembles 7 ensembles 









= 0.0005 


Ak — 


00025 




10^ 




1 0* 


± ^ 


10^ 


12^ 


K 














U.iObUU 






0.13(13) 


-0.16(22) 






0.15550 






0.01(10) 


-0.05(14) 






0.15600 


0.05(13) 


-0.26(32) 


-0.07(6) 


-0.06(16) 


0.13(8) 


0.08(13) 


0.15625 










0.08(6) 


0.01(14) 


0.15650 


-0.15(8) 


0.15(26) 


-0.01(5) 


-0.10(14) 


0.05(4) 


-0.02(13) 


0.15675 










0.02(3) 


0.06(10) 


0.15700 


-0.04(9) 


-0.16(16) 


-0.05(6) 


-0.01(9) 


0.04(2) 


-0.01(8) 


0.15725 










0.04(2) 


-0.09(5) 


0.15750 


0.11(7) 


0.05(5) 


0.00(4) 


0.10(11) 


0.04(2) 


-0.01(5) 


K 






m 








0.15500 


1.48(27) 


2.19(43) 


1.14(20) 


2.60(51) 






0.15550 






0.81(12) 


2.08(34) 






0.15600 


0.45(19) 


1.97(48) 


0.46(6) 


2.06(35) 


0.59(7) 


1.64(14) 


0.15625 










0.45(5) 


1.56(13) 


0.15650 


0.28(8) 


1.93(32) 


0.34(5) 


1.79(38) 


0.35(5) 


1.51(15) 


0.15675 










0.27(4) 


1.36(17) 


0.15700 


0.34(13) 


0.90(32) 


0.24(5) 


1.10(23) 


0.22(4) 


1.08(15) 


0.15725 










0.14(3) 


0.78(10) 


0.15750 


0.17(6) 


0.21(5) 


0.14(4) 


0.56(8) 


0.10(3) 


0.61(10) 
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Table 3: Swap acceptance data ai (3 — 5.6. 







6 ensembles 




7 ensembles 




0.155 <K< 0.1575 


0.156 <K< 0.1575 






Ak = 0.0005 




Ak = 0.00025 




8^ 


10^ 


12^ 


g4 20^ 12^ 


(A) 


0.63(1) 


0.41(1) 


0.22(1) 


0.82(1) 0.67(1) 0.54(1) 


lOOy/L^ 


0.98 


1.06 


0.99 


0.86 0.72 0.61 


lOy/L 


0.95 


1.06 


1.18 


0.69 0.72 0.73 
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standard HMC 



Tempered HMC 

6 ensembles 
0.155 <K< 0.1575 
Ak = 0.0005 



Tempered HMC 

7 ensembles 
0.156 <K< 0.1575 
Ak = 0.00025 




IIIi|IHI 



K = 0.156 



K = 0.15B 



Q t 





K = 0.156 



-5 



K = 0.1565 _ 



I L 



K = 0.1565 



K = 0.1565 





1 1 1 1 

K = 0.157 . 











5 - 

Q : 



K = 0.157 



K = 0.157 





K = 0.1575 _ 


5 

Q 




K = 0.1575 . 




^-^^A 













-5 







K = 0.1575 



rn 



traj. 



1000 



traj. 



1 000 



traj. 



Figure 1: Time series of Q for standard and tempered HMC on 10^ lattice at 
(3 = 5.6. The series for some intermediate K-values are not shown (see Table 
□ for full list). 
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standard HMC 



Tempered HMC 

6 ensembles 
0.155 <K< 0.1575 
Ak = 0.0005 



Tempered HMC 

7 ensembles 
0.156 <K< 0.1575 
Ak = 0.00025 





1 

K 


= 0.155 _ 








_ r 


' f ' ' '' '^' f 






1 





K = 0.156 



K = 0.15B 



K = 0.156 



K = 0.1565 _ 



i 



K = 0.1565 _ 



i4,i.!iiii'iii''H*i''^''i'iiiii',''' li'liii' 



' 



K = 0.1565 _ 



llllfl 



K = 0.157 



Q : 



K = 0.157 



Q t 



K = 0.157 



I "pBf,f,^i|iffl,^1^U(;11|lll!''^ 



K = 0.1575 



K = 0.1575 _ 
hi 



K = 0.1575 _ 



traj. 



1000 



traj. 



1 000 



traj. 



Figure 2: Time series of Q for standard and tempered HMC on 12^ lattice at 
(3 = 5.6. The series for some intermediate K-values are not shown (see Table 
□ for full list). 
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